function loglike_ts = loglike_model_cash_flow_timeseries_mv(param,data,data_zeroone,Tobs,nFirms,dt,Xhat1)
% Log-likelihood computed for each year t

%disp(param)

muP = param(1);
sigmaP = param(2);
alpha = param(3);
sigmaA = param(4);
rho  = param(5);

% TRANSITION  EQUATION: X_{t+1} = PhiX * X_tm1 + w_t,   w_t -> N(0,Q_t)
% MEASUREMENT EQUATION: Z_t = HZ * X_t + u_t,           u_t -> N(0,Omega_t)
PhiX = (1 + muP * dt); 
HZ = alpha * ones(nFirms,1) * dt; % nFirms x 1
Xhat = NaN(1,Tobs); % filtered state variable, 1 x 1
P = NaN(1,Tobs);
Q = NaN(1,Tobs); % variance of error in transition equation, transformed model
Xhat(1) = Xhat1; % true value of latent state at t=1
Z = data;
P(1) = sigmaP^2;

sigmaZZ = sigmaA * sqrt(dt); % stdev of measurement error 
sigmaXX = sigmaP * sqrt(dt); % stdev of transition equation error
rho_vec = rho*ones(nFirms,1);

Jvec = rho_vec(1) * sigmaZZ / sigmaXX * ones(nFirms,1);
HZstar = HZ + Jvec;
PhiXstar = - Jvec * PhiX;

loglike_ts = NaN(1,Tobs);
% number of firms can change each year

for t = 2 : Tobs % First observation is lost because of the one step ahead prediction in the Kalman filter
    % Kalman filter with non-zero correlation:
    Xhat_ttm1 = PhiX * Xhat(t-1); % prediction of X_t from t-1 to t
    Q(t) = sigmaXX^2 * Xhat(t-1)^2; % variance of transition error
	P_ttm1 = PhiX * P(t-1) * PhiX' + Q(t);
    
	% "Select matrix" of observed values: S_t is m_t x nFirms, and changes size at each t
    m_t = sum( data_zeroone(:,t) ); % number of firms/observations at date t
    S_t = zeros(m_t,nFirms); % initialize S_t, which is m_t x nFirms
    obs_t = find( data_zeroone(:,t) ); % find position(s) of observed values
    ind_t = sub2ind(size(S_t), [1:m_t]', obs_t);
    S_t(ind_t) = 1; % "Select matrix"

    HZstar_mv = S_t*HZstar; % adjust dimension
    PhiXstar_mv = S_t*PhiXstar; % adjust dimension    
    
    Z_ttm1 = HZstar_mv * Xhat_ttm1 + PhiXstar_mv * Xhat(t-1); % prediction of Z from t-1 to t
    Omegastar_mv = sigmaZZ^2 * Xhat(t-1)^2 * (1-rho_vec(1)^2) * eye(m_t); % cov-matrix (m_t x m_t) of measurement error
    
    F_ttm1 = HZstar_mv * P_ttm1 * HZstar_mv' + PhiXstar_mv * P(t-1) * PhiXstar_mv' + Omegastar_mv; % cov-matrix prediction error (Z_t - Z_ttm1)
    W = ( (P_ttm1 * HZstar_mv') / F_ttm1 )'; % Kalman gain, nFirms x 1

    Xhat(t) = Xhat_ttm1 + W' * ( S_t*Z(:,t) - Z_ttm1 ); % update of latent state variable X_t
    P(t) = P_ttm1 - W' * F_ttm1 * W;
    
    % Log-likelihood:
    logdetFttm1 = log( det(F_ttm1) );
    invFttm1 = eye(m_t) / F_ttm1; % = inv(F_ttm1);
    quadratic_form = ( S_t*Z(:,t) - Z_ttm1 )' * invFttm1 * ( S_t*Z(:,t) - Z_ttm1 );
    ntlog2pi = m_t * log(2*pi); % Note: number of firms changes at each year t
    loglike_ts(t) = - 0.5 * ( ntlog2pi + logdetFttm1 + quadratic_form ); 
    
end

loglike_ts = loglike_ts(2:end); % first observation is lost



